
/***
 * Funtion for time-density curve measurement in ROIs on DSAs.
 * With the optional use of vessel masking.
 ***/
 
#include <iostream>
#include <vector>

using namespace std;

#include "itkImageRegionConstIterator.h"

#include "dsaio.h"
#include "roi_mask.h"


/**
 * Performes densitometry.
 * @return vector<double> average density inside polygonal ROI on each frame
 */
vector< double > dmetry ( char *dsa_file, char *roi_file )
{ 
	// read ROI points
	ifstream input_roi;
	input_roi.open ( roi_file );
	vector< int > roi_x, roi_y;
	int n, x, y;
	while ( input_roi >> n )
	{
		input_roi >> x;
		input_roi >> y;
		if ( ! input_roi.bad() )
		{
			roi_x.push_back( x );
			roi_y.push_back( y );
		}
	}
	input_roi.close();

	// compute ROI mask
	vector< unsigned char > mask = roi_mask ( roi_x, roi_y );

	// determine bounding box of ROI
	int min_x = 20000;
	int max_x = -1;
	int min_y = 20000;
	int max_y = -1;
	for ( int i = 0; i < roi_x.size(); ++ i )
	{
		if ( roi_x[ i ] < min_x ) min_x = roi_x[ i ];
		if ( roi_x[ i ] > max_x ) max_x = roi_x[ i ];
		if ( roi_y[ i ] < min_y ) min_y = roi_y[ i ];
		if ( roi_y[ i ] > max_y ) max_y = roi_y[ i ];
	}
	int width = max_x - min_x;
	int height = max_y - min_y;

	// cerr << "densitometry" << endl;
	// densitometry
	vector< double > dcurve;
	typedef itk::ImageRegionConstIterator< DsaFrame > DsaConstIterator;
	int N = dsa_frame_count ( dsa_file );

	for ( n = 0; n < N; ++ n )// iterate through all frames
	{
		DsaFrame::Pointer frame = dsa_get_frame ( dsa_file, n );
		DsaFrame::RegionType roi;
		DsaFrame::RegionType::IndexType roiStart;
		DsaFrame::RegionType::SizeType roiSize;
		roiStart[ 0 ] = min_x;
		roiStart[ 1 ] = min_y;
		roiSize[ 0 ] = width;
		roiSize[ 1 ] = height;
		roi.SetSize( roiSize );
		roi.SetIndex( roiStart );
		DsaConstIterator it ( frame, roi );
		int i = 0;
		double avg = 0.0;
		double denom = 0.0;
		for ( it.GoToBegin(); ! it.IsAtEnd(); ++ it )// iterate pixels
		{
			if ( (int)mask[ i ] != 0 )
			{
				avg += (double)( it.Get() );
				denom += 1.0;
			}
			++ i;
		}
		avg /= denom;

		dcurve.push_back ( avg );
	}

	return dcurve;
}

/**
 * Uses a vessel mask image file (same size as DSA), where black (=0) points
 * indicate measurement points. All positions, not zero on vessel mask will
 * be ignored during densitometry.
 */
vector< double > dmetry ( char *dsa_file, char *roi_file, char *vm_file )
{ 
	// read ROI points
	ifstream input_roi;
	input_roi.open ( roi_file );
	vector< int > roi_x, roi_y;
	int n, x, y;
	while ( input_roi >> n )
	{
		input_roi >> x;
		input_roi >> y;
		if ( ! input_roi.bad() )
		{
			roi_x.push_back( x );
			roi_y.push_back( y );
		}
	}
	input_roi.close();

	// compute ROI mask
	vector< unsigned char > mask = roi_mask ( roi_x, roi_y );

	// determine bounding box of ROI
	int min_x = 20000;
	int max_x = -1;
	int min_y = 20000;
	int max_y = -1;
	for ( int i = 0; i < roi_x.size(); ++ i )
	{
		if ( roi_x[ i ] < min_x ) min_x = roi_x[ i ];
		if ( roi_x[ i ] > max_x ) max_x = roi_x[ i ];
		if ( roi_y[ i ] < min_y ) min_y = roi_y[ i ];
		if ( roi_y[ i ] > max_y ) max_y = roi_y[ i ];
	}
	int width = max_x - min_x;
	int height = max_y - min_y;

	// cerr << "densitometry" << endl;
	// densitometry
	vector< double > dcurve;
	typedef itk::ImageRegionConstIterator< DsaFrame > DsaConstIterator;
	int N = dsa_frame_count ( dsa_file );

	for ( n = 0; n < N; ++ n )// iterate through all frames
	{
		DsaFrame::Pointer frame = dsa_get_frame ( dsa_file, n );
		DsaFrame::Pointer vm = dsa_get_frame ( vm_file, n );
		DsaFrame::RegionType roi;
		DsaFrame::RegionType::IndexType roiStart;
		DsaFrame::RegionType::SizeType roiSize;
		roiStart[ 0 ] = min_x;
		roiStart[ 1 ] = min_y;
		roiSize[ 0 ] = width;
		roiSize[ 1 ] = height;
		roi.SetSize( roiSize );
		roi.SetIndex( roiStart );
		DsaConstIterator it ( frame, roi );
		DsaConstIterator vmit ( vm, roi );
		int i = 0;
		double avg = 0.0;
		double denom = 0.0;
		for ( it.GoToBegin(), vmit.GoToBegin(); ! it.IsAtEnd(); ++ it, ++ vmit )
		{
			if ( (int)mask[ i ] != 0 && vmit.Get() == 0 )
			{
				avg += (double)( it.Get() );
				denom += 1.0;
			}
			++ i;
		}
		avg /= denom;

		dcurve.push_back ( avg );
	}

	return dcurve;
}
